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Abstract 

A method for the combination of microscopic and macroscopic simulations is devel- 
oped which is based on the invariance of the macroscopic relative to the microscopic 
dynamics. The method recognizes the onset and breakdown of the macroscopic de- 
scription during the integration. We apply this method to the case of ferrofluid 
dynamics, where it switches between direct Brownian dynamics simulations and 
integration of the constitutive equation. 
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1 Introduction 



The understanding of the macroscopic dynamics from the underlying micro- 
scopic time evolution is a central issue of non-equilibrium statistical mechan- 
ics. The massive use of computer simulations over the last years has led to 
new approaches to this very old problem. Among others, we mention Legendre 
integrators [1,2,3], the CONNFESSIT method [4], adaptive mesh refinement 
and multiscale modeling [5,6]. The last two methods do not require knowledge 
of the macroscopic equations. On the other hand, there has been much effort 
to derive approximate macroscopic equations from the microscopic dynamics, 
which yield reliable results under various circumstances (see e.g. Ref. [7] for 
an overview of constitutive equations for polymer liquids). Here, we follow 
the approach proposed in Ref. [8] to combine microscopic and macroscopic 
simulations in a combined integration scheme which recognizes the onset and 
breakdown of the chosen macroscopic description during the simulation. Note, 
that the breakdown of a chosen macroscopic description does not imply a sim- 
ilar breakdown of other, improved macroscopic descriptions. Instead of im- 
proving the macroscopic equations, which is the aim of many works on closure 
approximations (see e.g. [1,2] and references therein), we here keep the chosen 
macroscopic description and use it as long and as frequently in the simulation 
as possible. While this integration scheme was used in Ref. [8] to detect the 
onset of the macroscopic description, we here present the full scheme that 
switches back and forth between microscopic and macroscopic simulations. 
We apply this scheme to well-known models of ferrofluid dynamics where it 
decides between direct Brownian dynamics simulations and integration of the 
constitutive equation. 



2 Invariance principle and combined integration scheme 

In order to keep the paper self-contained, we briefly summarize the main 
ideas of the combined integration scheme based on the invariance principle 
proposed in Refs. [8,9]. We assume a given microscopic description of the 
system, where the microscopic variables are denoted by /. The microscopic 
dynamics is specified by the vector field J, 
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In addition, we assume that the set of macroscopic variables M has been 
chosen. Typically, f(x) is the distribution function over the set of microscopic 
coordinates x and M contains low-order moments of /. In this case, the macro- 
scopic variables are linear functionals of the microscopic distribution function, 
M = Jdxm(x)f(x). Although the method can be applied to more general 
situations, we here limit ourselves to this case for the sake of clarity. 

The reduced or macroscopic description assumes not only closed-form macro- 
scopic equations M(M), but also a family of canonical distribution functions 
/m(x) [1,2]. The canonical distribution functions satisfy the consistency rela- 
tion M = jdxm(x)f-M.(x). Then, the macroscopic dynamics is given by 

M(M) = jdxm{x)J{f m {x)) (2) 



Different routes to the construction of /m have been proposed. In many ap- 
plications, the dynamic system Eq. (1) is equipped with a Lyapunov function 
S (the entropy, free energy, etc.), and the canonical distribution functions /m 
are conditional maximizers of S subject to fixed M [8]. 

In order to estimate the accuracy of the macroscopic description we define the 
defect of invariance A M (x) as the difference of the microscopic and macro- 
scopic time derivative, 

AmW = J(/m)-^-M(M). (3) 



By construction, fdxm(x)A^/[(x) = 0. If the defect of invariance Am(s) van- 
ishes for all admissible values of M, then the reduced description is called 
invariant and the family /m represents the invariant manifold in the space 
of the microscopic variables. The invariant manifold is relevant if it is stable. 
Exact invariant manifolds are known only in very few cases. Corrections to 
the manifold /m through minimizing Am is part of the so-called method of 
invariant manifolds [10,11]. 

Here, we exploit the invariance principle in a different way. Let f(x; t\to) denote 
the microscopic variables at time t for given initial conditions at time to < t. 
The values of the macroscopic variables at time t are given by M(t|t ) = 
Jdxm(x)f(x;t\t ). On the other hand, the solution of the macroscopic equa- 
tions (2) with corresponding initial conditions gives M*(i|£o)- We denote with 
|| A || the value of the defect of invariance (3) with respect to some norm || • || 
and e > a fixed threshold value. If at time t the defect of invariance satisfies 

||A M (t|t )|| < e, (4) 
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it is said that the macroscopic description sets on, since the reduced description 
is sufficiently accurate. However, if 

||A M *(t|t )|| > e, (5) 



the macroscopic description breaks down since the accuracy of the macroscopic 
dynamics is insufficient. Therefore, the evaluation of the defect of invariance 

(3) on the current solution either to the macroscopic or to the microscopic 
dynamics and checking Eqs. (4) and (5) we can decide whether integration of 
the macroscopic dynamics is sufficiently accurate or not. 

This information is used in the combined integration scheme to switch be- 
tween microscopic and macroscopic simulations. The scheme is sketched in 
Fig. 1. Suppose at time £ the microscopic dynamics is integrated for given 
initial condition. The integration is continued until at time t\ the inequality 

(4) is satisfied. At this point, the macroscopic dynamics is started with the 
actual values of the macroscopic variables, M*(£i) = M(£i|£ ). The macro- 
scopic dynamics is integrated until the macroscopic description breaks down 
at a later time £2, which is signaled by ||AM*(t 2 |ti)|| > e - At this time it is 
necessary to switch back from the macroscopic to the microscopic simulations 
in order to achieve the required accuracy. The initial condition for the mi- 
croscopic simulation at time £2 is obtained from the macroscopic description, 
f(x, £2) = f-M(t 2 )( x )- Then, the microscopic dynamics is integrated until the 
macroscopic description sets on etc. In the sequel, we demonstrate this scheme 
for the case of ferrofluid dynamics. 



3 Kinetic models of ferrofluid dynamics 



Ferrofluids are stable suspensions of nano-sized ferromagnetic colloidal parti- 
cles in a suitable carrier liquid [12]. These fluids attract considerable interest 
due to their peculiar behavior, such as the magnetoviscous effect, the depen- 
dence of the viscosity coefficients on the magnetic field [13]. 

We here consider the kinetic model of ferrofluid dynamics proposed in Refs. [14,15]. 
In this model, the ferromagnetic particles are assumed to be identical, mag- 
netically hard ferromagnetic monodomain particles. It is further assumed that 
the particles are of an ellipsoidal shape with axes ratio r and that the magnetic 
moment \i is oriented parallel to the symmetry axes of the particle. Let /(u) 
denote the orientational distribution function to find a ferromagnetic parti- 
cle with the orientation u, where u is a vector on the three-dimensional unit 
sphere. In the general notation of Sec. 2, the microscopic coordinates x are the 
orientations u and f(x) is the orientational distribution function /(u). The 
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normalized macroscopic magnetization M is given by M = / d 2 u u/(u), where 
/ d 2 u denotes integration over the three-dimensional unit sphere. The normal- 
ization is performed with the saturation magnetization M sat = n/x, where n 
denotes the number density. In the presence of a local magnetic field H and a 
velocity field v, the dynamics is given by [12,16] 

d 

—f = -C ■ {(Q x u + Bu x D ■ u - D r h x u)/} + D r C ■ Cf (6) 



In Eq. (6) we have introduced the rotational operator C = u x d/du, the 
vorticity Q = V x v/2, the symmetric velocity gradient 2D = Vv + (Vv) T , 
and the so-called shape factor B = (r 2 — l)/(r 2 + 1). The rotational diffu- 
sion coefficient D r defines the rotational relaxation time r = (2D r )~ l . The 
dimensionless magnetic field is defined by h = /xH/ZcbT, where k-g and T 
denote Boltzmann's constant and temperature, respectively. The equilibrium 
distribution 

/h(u) = - k exp [h • u] (7) 
47r smh(/i) 



is the stationary solution to the kinetic equation (6) in the absence of flow. 
From Eq. (7), the equilibrium magnetization M cq is found to be given by 
M cq /M sat = Li(h)h/h, where h = Vh • h is the Langevin parameter and 
Li(x) = coth(x) — x' 1 denotes the Langevin function. 

Except for special cases, exact solutions to the kinetic equation (6) are un- 
known and closed form equations for the magnetization cannot be derived 
exactly from Eq. (6). In order to solve the closure problem, the authors of 
Ref. [17] have suggested to use the family of equilibrium distributions (7), 
where the magnetic field h is replaced by an effective field Thus, the non- 
equilibrium magnetization is given by M/M sat = Li(£)n, where we have intro- 
duced the norm £ of the effective field £ and n = This so-called Effective 
Field Approximation (EFA) is a particular instance of the quasi-equilibrium 
or maximum entropy approximation. It is derived from extremizing the en- 
tropy functional S[f] = — Jd 2 u /(u) ln[/(u)// h (u)] subject to the constraints 
of fixed normalization and fixed values of magnetization M. Therefore, the 
set of macroscopic variables contains only the magnetization M in the present 
case. For more details on the use of quasi-equilibrium approximation in the 
context of complex fluids see e.g. [1,2,3]. The macroscopic equation (2) be- 
comes the magnetization equation 



M = fl xM + B 



5 



-2D r M + D r 1 - 




h- 



^(0 



MO 



(h- n)M 



(8) 



Functions Li(x) are defined recursively by L i+ i(x) = Li-i(x) — (2i + l)Li(x) / x, 
with L (x) = 1 and Li(x) the Langevin function. The accuracy of the EFA 
has been discussed, e.g., in Refs. [17,18,19,20]. 

The defect (3) of this approximation can be calculated explicitly. Note, that 
/ d 2 u A M as well as / d 2 u uA M vanish identically by construction. Therefore, 
we estimate the accuracy of the EFA through information about the dynamics 
of the next higher order moment that is not included in the macroscopic 
description and consider the matrix Am = / d 2 u uuAm- As a suitable norm we 



use the matrix norm || A|| = y Sq^Aq^) 2 . The matrix Am can be represented 

by 

A M = dil + d 2 nn + d 3 (hn + nh) + d 4 ~D + d 5 (D • nn + nn ■ D), (9) 

where the coefficients di are defined in the Appendix A. 

4 Numerical implementation 

The microscopic dynamics (6) is integrated by Brownian dynamics simulations 
of the corresponding Ito stochastic differential equation 



The projector perpendicular to U t is denoted by P t = (1 — U t U t ) and W t is a 
three-dimensional Wiener process [21]. Using Ito's formula, it is verified that 
Eq. (10) conserves the normalization of Ut- Eq. (10) is integrated numerically 
by a weak first-order scheme that guarantees the normalization of the random 
unit vector U t [21]. An ensemble of 10 5 random vectors U t is used in the 
simulation in order to ensure accurate ensemble averages. A constant time 
step At/r = 10~ 3 is used throughout. 

The macroscopic equation (8) is integrated directly in terms of the macroscopic 
variables M. The effective field £ is calculated by £ = Lj" 1 (M), where M 
denotes the norm of M and Lf 1 (a;) denotes the inverse Langevin function. 
The latter is evaluated numerically by the Newton-Raphson Method. In order 
to treat microscopic and macroscopic dynamics approximately on the same 
footing, we used a first-order explicit Euler scheme with the same time step 
At/r = 10~ 3 to integrate Eq. (8) numerically. More accurate schemes for the 




dU t =P t ■ [(ft xU t + BB-U t + h)dt + dW t ] - U t 
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macroscopic equation, such as Runge-Kutta methods, could be employed as 
well. 

At every time step of the microscopic or macroscopic integration, the norm of 
the defect of invariance || A|| is evaluated from Eq. (9). If in the course of the 
microscopic integration the inequality ||A|| < e is satisfied at time i, then the 
microscopic simulation is stopped and the macroscopic integration of Eq. (8) 
is started with the initial condition M(i) calculated from the microscopic 
ensemble average. On the other hand, if at time i during the integration of 
the macroscopic equation ||A|| > e is fulfilled, the microscopic simulation is 
started with initial condition /(u;i) = /m(*)(u). Here, we use the rejection 
method [22] in order to generate an ensemble of random unit vectors Ut, 
which are distributed according to /M(t)( u )- 



5 Results 

The combined integration scheme described above has been implemented and 
run for a number of different values of magnetic field and velocity gradients. 
We limited ourselves to the case of plane shear flow. Previous investigations 
showed that the EFA provides a good approximation to the kinetic equation 
(6) for moderate values of axes ratios r [20]. We here choose a value of r = 5 
in the sequel, where the accuracy of the EFA is not as good as for smaller 
values. 

First, we consider the dynamics in the absence of any velocity gradients where 
a constant magnetic field h is applied during the time interval t\ < t < if. 
Comparison of the BD simulation results and the EFA for h = 2 and h = 5 
with i; = r, if = 4r is shown in Fig. 2. From Fig. 2 we observe that the EFA 
provides a very good approximation for this case. Deviations of the EFA from 
the results of the BD simulation are shown in Fig. 3 together with the values 
of the norm of the defect of invariance. As can be seen from Fig. 3, the defect 
is very sensitive to deviations of the EFA from the BD results. Upon closer 
inspection one recognizes from Fig. 3 that the maxima of the defect preceed 
the extrema of the deviations, reflecting the fact that the defect of invariance 
is sensitive to deviations in the time derivative. Thus, the defect signals the 
inaccuracy of the time derivative of the EFA and allows for a switch to the 
BD simulation before the values of the magnetization become inaccurate. The 
result of the combined integration scheme with e = 0.2 is shown in Fig. 4. 
Within the boxed regions, the norm of the defect of invariance ||A|| > e and 
the BD simulation is performed while otherwise the EFA is integrated. In 
the inset of Fig. 4, the result of the EFA, Fig. 2, is shown for comparison. 
Clearly, the combined integration improves the comparison of the EFA to the 
BD simulation. 
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Next, we consider the magnetization dynamics in the presence of flow. Starting 
with an isotropic initial distribution in the absence of magnetic fields and 
velocity gradients, a constant magnetic field h is applied during the interval 
ti < t < if, while the magnetic field is absent outside this interval. In addition, 
a plane Couette flow with velocity field v = (2^y, 0, 0) with a constant shear 
rate 7 is applied during the interval t\ < t < t 2 while no shear is applied for 
t < t\ and t > t 2 . Fig. 5 shows the magnetization dynamics for h = 5, T7 = 2, 
ti — t, tf — 6r, t x = 3r, and t 2 = St. One observes from Fig. 5 that the EFA 
predicts qualitativly correct behavior but fails to give accurate results as long 
as the shear flow is applied. Fig. 6 shows the deviation of the EFA from the BD 
simulation results together with the norm of the defect of invariance. Except 
for initial transient dynamics when the magnetic field is suddenly applied (see 
Fig. 3) the EFA is accurate in the absence of shear flow and ||A|| is small. 
During application of the shear flow, however, the predictions of the EFA 
are less accurate as is monitored by ||A||. In Fig. 7 we show the result of 
the combined integration scheme for the present situation where e = 0.2. The 
agreement between the combined scheme with the BD simulation is very good. 
Due to the low accuracy of the EFA, the macroscopic equation (8) is integrated 
only in the boxed regions while otherwise BD simulations are performed. 

Finally, we consider the magnetization dynamics in a constant magnetic field 
h and a plane Couette flow v = (2 A /y, 0, 0) with oscillatory shear rate 7 = 
wfo cos(cjt). The magnetic field is oriented in the flow direction. Equilibrium 
initial conditions are chosen. Due to the flow, a nonequilibrium magnetization 
component M y arises which oscillates with frequency u, while the magnetiza- 
tion component in the magnetic field direction M x oscillates with frequency 
Iuj. From Fig. 8 we observe that the predictions of the EFA are less reliable 
compared to the situation without flow. The deviation of both magnetiza- 
tion components are shown in Fig. 9 together with the defect of invariance. 
The deviation of the magnetization components of the macroscopic dynamics 
(EFA) from the microscopic dynamics (BD) are seen to oscillate as well. While 
the deviation of M y oscillates around zero, the macroscopic dynamics always 
overpredicts the values of M x in the present case. The oscillations in the de- 
viation of M x and M y are reflected in the oscillations of ||A||. As before, the 
maxima of ||A|| occur before the maximum deviations of the magnetization 
are observed. Note that ||A|| always exceeds a certain value e , which signals 
the inaccuracy of the EFA in the present case. Fig. 10 shows the result of 
the combined integration scheme with e = 0.2 for the same conditions as in 
Fig. 8. In the boxed regions ||A|| < e and the macroscopic (EFA) dynamics 
is integrated while otherwise BD simulations are performed. The agreement 
between combined integration and BD simulation is very good. However, due 
to the limited accuracy of the EFA for this case only a limited fraction of 
the total integration time is preformed with the EFA. Fig. 11 shows the same 
situation as in Fig. 8, but the shear flow was stopped at time t = 6.5r while 
the magnetic field remained unchanged. Without the terms arising from the 
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velocity gradient, || A|| drops below the threshold value e and, within the com- 
bined scheme, the relaxational dynamics for t > 6.5r is integrated by the EFA. 
Good agreement with the result of full BD simulation is found. 



6 Conclusions 

The present approach of combined microscopic and macroscopic simulation ex- 
ploits the invariance of the microscopic relative to the macroscopic dynamics. 
It is applicable whenever a macroscopic description to an underlying micro- 
scopic dynamics is given. While previous studies considered the case of dilute 
polymer solutions [8,9], the combined integration scheme is illustrated here 
for the case of ferrofluid dynamics. 

The microscopic dynamics is integrated only if the defect of invariance ex- 
ceeds a certain threshold value e. The full microscopic simulation is recovered 
for e = while e — > oo corresponds to the macroscopic dynamics. Thus, 
the combined integration improves the accuracy of the macroscopic descrip- 
tion, where the improvement depends on the choice of e. At the same time, 
the combined integration saves CPU time since the macroscopic simulation 
is employed whenever possible. The amount of CPU time that can be saved 
by the combined integration scheme for a given value of e depends on the 
quality of the macroscopic description for the corresponding situation. For 
the conditions of Fig. 8, a sample study of the relative error R x of the mag- 
netization M x as a function of elapsed CPU time is shown in Fig. 12. The 
relative error with respect to the result of the BD simulation is defined as 
R\ = J- J2f l [(M x (tj) — M® D (tj))/ M® D (tj)] 2 , where N t denotes the total num- 
ber of integration time steps. For a better comparison, all data shown in Fig. 12 
are obtained with the same PC with a P4 processor. From Fig. 12 we observe 
that the relative error decreases with decreasing e while the time the micro- 
scopic simulation is integrated in the combined scheme increases and thus the 
required CPU time increases. Overall, we observe that the relative error de- 
creases almost linearly with elapsed CPU time. Note, that R x = does not 
correspond to the exact result but to the BD simulation. 
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A Coefficients of defect of invariance 



The coefficients di in Eq. (9) contain contributions from Brownian motion, the 
magnetic field and the symmetric velocity gradient. In particular, 

* = XMO + cm (i~)-b (2^ + (D: nn)(A ' 1) 

d 2 = -6(L 2 ({) + c(0) (l - ^) + 2 UiiOljJgj " j] " MO) n ■ h 

+ IH«- 4 ^l +9 ^l) D:m (A ' 2) 



(A.3) 
(A.4) 
(A.5) 



where 



c(0 = ^[£ 2 (0 - MO 2 ] (A.6) 



and the total derivative of the Langevin function can be expressed by 



10 



References 



[I] P. Ilg, I. V. Karlin, H. C. Ottinger, Canonical distribution functions in polymer 
dynamics: I. Dilute solutions of flexible polymers, Physica A 315 (2002) 367-385. 

[2] P. Ilg, M. Kroger, I. V. Karlin, H. C. Ottinger, Canonical distribution functions 
in polymer dynamics: II. liquid-crystalline polymers, Physica A 319C (2003) 
134-150. 

[3] A. N. Gorban, P. A. Gorban, I. V. Karlin, Legendre integrators, post-processing 
and quasiequilibrium, J. Non-Newton. Fluid Mech., this issue. 

[4] M. Laso and H. C. Ottinger, Calculation of viscoelastic flow using molecular 
models - the CONNFESSIT approach, J. Non-Newton. Fluid Mech. 47 (1993) 
1-20. 

[5] R. M. Jendrejack, J. J. de Pablo, M. D. Graham, A method for multiscale 
simulation of flowing complex fluids, J. Non-Newton. Fluid Mech. 108 (2002) 
123-142. 

[6] C. I. Siettos, M. D. Graham, I. G. Kevrekidis, Coarse Brownian dynamics 
for nematic liquid crystals: Bifurcation, projective integration, and control via 
stochastic simulation, J. Chem. Phys. 118 (2003) 10149-10156. 

[7] R. B. Bird, J. M. Wiest, Constitutive equations for polymeric liquids, Ann. Rev. 
Fluid Mech. 27 (1995) 169-193. 

[8] A. N. Gorban, I. V. Karlin, P. Ilg, H. C. Ottinger, Corrections and enhancements 
of quasi-equilibrium states, J. Non-Newton. Fluid Mech. 96 (2001) 203-219. 

[9] I. V. Karlin, P. Ilg, H. C. Ottinger, Invariance principle to decide between 
micro and macro computations in Recent Developments in Mathematical and 
Experimental Physics, Volume C: Hydrodynamics and Dynamical Systems. 
p. 43-50, F. Uribe (Ed.), Kluwer, Dordrecht, 2002. 

[10] A. N. Gorban and I. V. Karlin, Thermodynamic parameterization, Physica A 
190 (1992) 393-404. 

[II] A. N. Gorban, I. V. Karlin, A. Y. Zinovyev, Constructive methods 
of invariant manifolds for kinetic problems, preprint 2003, available at 
http://www.ihes.fr/PREPRINTS/M03/Resu/resu-M03-50lrtrnT 

[12] E. Blums, A. Cebers, M. M. Maiorov, Magnetic Fluids, de Gruyter, Berlin, 
1997. 

[13] M. Kroger, P. Ilg, S. Hess, Magnetoviscous model fluids, J. Phys. 
Condens. Matter 15 (2003) S1401-S1423. 

[14] M. A. Martsenyuk, Viscosity of a suspension of ellipsoidal ferromagnetic 
particles in a magnetic field, J. Appl. Mech. Tech. Phys. 14 (1973) 564-566. 



11 



[15] A. Cebers, Simulation of the magnetic rheology of a dilute suspension of 
ellipsoidal particles in a numerical experiment, Magnetohydro dynamics 20 
(1984) 349-354. 

[16] S. Hess, Fokker-Planck-Equation approach to flow alignment in liquid crystals 
Z. Naturforsch. 31a (1976) 1034. 

[17] M. A. Martsenyuk, Yu. L. Raikher, M. I. Shliomis, On the kinetics of 
magnetization of suspension of ferromagnetic particles, Sov. Phys. JETP 38 
(1974) 413-416. 

[18] P. Ilg, M. Kroger, S. Hess, Orientational Order 
Parameters and Magnetoviscosity of Dilute Ferrofluids, J. Chem. Phys. 116 
(2001) 9078-9088. 

[19] P. Ilg and M. Kroger, Magnetization dynamics, rheology, and an effective 
description of ferromagnetic units in dilute suspension, Phys. Rev. E 66 (2002) 
021501. 

[20] P. Ilg, M. Kroger, S. Hess, A. Yu. Zubarev, Dynamics of colloidal suspensions 
of ferromagnetic particles in plane Couette flow: comparison of approximate 
solutions with Brownian dynamics simulations, Phys. Rev. E 67 (2003) 061401. 

[21] H. C. Ottinger, Stochastic Processes in Polymeric Fluids. Springer, Berlin, 1996. 

[22] W. H. Press, S. A. Teukolsky, W. T. Vellering, B. P. Flannery, Numerical Recipes 
in Fortran. Cambidge University Press, 2nd ed. 1992. 



12 



Figure captions: 



Fig. 1 Sketch of the combined integration scheme. The macroscopic dynamics 
is integrated whenever the norm of the defect of invariance ||A|| is smaller 
than some fixed threshold value e. Otherwise, the microscopic dynamics is 
integrated. 

Fig. 2 Magnetization dynamics as a function of reduced time t/r in the 
absence of velocity gradients. A constant magnetic field h was applied during 
the time interval 1 < t/r < 4, while the magnetic field was switched off 
outside this interval. Circles and squares are the results of the BD simulation 
for h = 5 and h = 2, respectively, while solid and dashed line are the 
corresponding predictions of the EFA. The inset shows the comparison for 
h = 5 on a finer scale. 

Fig. 3 Deviation of normalized magnetization M x /M sait calculated from BD 
simulation and the EFA for h = 5 (circles) and h = 2 (squares) for the same 
conditions as in Fig. 2. Solid and dashed lines are the defect of invariance as 
calculated from the matrix norm of Eq. (9) for h = 5 and h = 2, respectively. 
For better visibility, the matrix norm was multiplied by a factor 0.1. 

Fig. 4 Magnetization dynamics as a function of reduced time t/r for the 
same condition as in Fig. 2. Circles and squares are the result of the BD 
simulation for h — 5 and h = 2, respectively, while solid and dashed lines 
are the results of the combined integration scheme with e = 0.2 for h = 5 
and h = 2, respectively. Within the boxed regions (indicated by the shading 
in the upper part), ||A|| > e and the BD simulation is performed, otherwise 
the EFA is integrated. The inset shows the comparison for h = 5 on a finer 
scale, where the dashed-dotted line is the result of the EFA. 

Fig. 5 Magnetization dynamics as a function of reduced time t/r in a con- 
stant magnetic field h = 5 for 1 < t/r < 6 and steady shear flow with shear 
rate T7 = 2 for 3 < t/r < 8, where the magnetic field is oriented in the 
gradient direction. No magnetic field and no shear flow is applied outside 
the mentioned time intervals. Circles and squares represent the results of 
the BD simulation for M x /M sat and M y /M sat , respectively, while solid and 
dashed lines are the corresponding result of the EFA. 

Fig. 6 Deviation of normalized magnetization M x /M sat (circles) and M y /M sat 
(squares) calculated from BD simulation and the EFA for the same condi- 
tions as in Fig. 5. The solid line is the defect of invariance as calculated 
from the matrix norm of Eq. (9). For better visibility, the matrix norm was 
multiplied by a factor 0.1. 

Fig. 7 Magnetization dynamics as a function of reduced time t/r for the same 
conditions as in Fig. 5. Circles and squares represent the result of full BD 
simulation, the dashed-dotted line corresponds to the EFA and full lines are 
the result of the combined integration scheme, where the EFA is integrated 
within the boxed regions while otherwise BD simulations are performed. 

Fig. 8 Magnetization dynamics as a function of reduced time t/r for in- 
ception of oscillatory shear flow with frequency tuj = 2 and amplitude 
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7o = 0.5. The magnetic field is oriented in flow direction with h = 2. Circles 
and squares represent the results of the BD simulation for M x /M SSut and 
M y /M sat , respectively, while solid and dashed lines are the corresponding 
result of the EFA. 

Fig. 9 Deviation of normalized magnetization M x /M &&t (circles) and M y /M ssX 
(squares) calculated by BD simulation and from EFA as functions of time 
t/r. Also shown is the norm of defect of invariance ||A|| (solid line), which 
is multiplied by a factor 0.1 for better visibility. The same flow conditions 
as in Fig. 8 are considered. 

Fig. 10 Magnetization dynamics as a function of reduced time t/r for the 
same conditions as in Fig. 8. Symbols are the result of the BD simulation. 
Solid and dashed lines correspond to the combined integration, where the 
EFA is integrated within the boxed regions (indicated by the shading in the 
upper part) and the BD simulation is preformed outside. 

Fig. 11 Magnetization dynamics as a function of reduced time t/r for the 
same conditions as in Fig. 8, but where the shear flow was stopped at time 
t = 6.5t. Symbols are the result of the BD simulation. Solid and dashed 
lines correspond to the combined integration, where the EFA is integrated 
within the boxed regions (indicated by the shading in the upper part) and 
the BD simulation is preformed outside. Dashed-dotted lines are the result 
of the EFA. 

Fig. 12 Relative error R x defined in the text as a function of CPU time in 
seconds on a logarithmic scale. The same conditions as in Fig. 8 are con- 
sidered. The number above the filled symbols are the corresponding values 
of e. Solid lines are guides to the eye. Different values of e decreasing from 
e = 1 (EFA) to e = (BD simulation) have been chosen in the combined 
integration scheme in order to obtain increasingly more accurate results for 
M x . 



14 



macro 



to 



t 2 



micro 



l|A||>e 



||A||<6 ||A||>e 



Fig. 1. 



15 




reduced time 



Fig. 2. 



16 




Fig. 3. 



17 




Fig. 4. 



18 




20 




Fig. 7. 



21 



0.1 




Fig. 9. 



23 



3 




reduced time 



Fig. 10. 



24 



to 





6000 



CPU time [sec] 
Fig. 12. 



26 



